
install.packages("tseries"); install.packages("egcm");install.packages("car")
library(tseries); library(egcm); library(car)


# setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Main_Text/Figure_1_&_2")



data_q<-read.csv(file="data_fig_1_2.csv")


############## Figure 1 ##################


mean_ts<-tapply(data_q$Area,data_q$Year,mean,na.rm=T)

median_ts<-tapply(data_q$Area,data_q$Year,median)

quantile_25ts<-tapply(data_q$Area,data_q$Year,quantile,probs=c(.25))

quantile_75ts<-tapply(data_q$Area,data_q$Year,quantile,probs=c(.75))




mean_Lts<-tapply(data_q$logArea,data_q$Year,mean,na.rm=T)

median_Lts<-tapply(data_q$logArea,data_q$Year,median)

quantile_25Lts<-tapply(data_q$logArea,data_q$Year,quantile,probs=c(.25))

quantile_75Lts<-tapply(data_q$logArea,data_q$Year,quantile,probs=c(.75))




#### Top Panel ####

years<-as.numeric(as.character(unique(data_q$Year)))



pdf("Figure_1_Top_Panel.pdf",width=7, height=5)


plot(years, mean_ts,col="black",lty=1,type="n",,ylim=c(0,75000),ylab="Square Kilometers",main="Untransformed Data",xlab="",axes=F)

grid()

lines(years,median_ts,xlab="Year",lwd=2)

lines(years, mean_ts,col="dark grey",lty=1,lwd=2)  

lines(years, quantile_25ts,lty=2,lwd=1)  

lines(years, quantile_75ts,lty=2,lwd=1)  





text(1205,24500,"75th",cex=.5)
text(1090,0,"25th",cex=.5)
text(1790,57500,"Mean",cex=1)
text(1110,11000,"Median",cex=1)

axis(2,seq(0,80000,by=20000))
axis(1)


dev.off()



#### Lower Panel ####

pdf("Figure_1_Lower_Panel.pdf",width=7, height=5)


plot(years,median_Lts,type="n",ylim=c(4,11),xlab="Year",ylab="log(Square Kilometers)",lwd=2,main="Transformed Data",axes=F,width=7, height=5)

grid()

lines(years,median_Lts,xlab="Year",lwd=2)

lines(years, mean_Lts,col="dark grey",lty=1,lwd=2)  

lines(years, quantile_25Lts,lty=2,lwd=1)  

lines(years, quantile_75Lts,lty=2,lwd=1)  




text(1350,9.5,"75th",cex=.5)

text(1500,4.5,"25th",cex=.5)

text(1790,6.4,"Mean",cex=1)

text(1790,5.55,"Median",cex=1)


axis(1)
axis(2)


dev.off()


######################### Dickey Fuller Test Footnote 39 ####################



out<-lm(mean_Lts~median_Lts-1)


dickey_fuller<-adf.test(resid(out),alternative="stationary",k=0)


print("dickey fuller test footnote 39")

print(dickey_fuller)


############## Figure 2 ##################



pdf(file="Figure_2.pdf")

par(mfrow=c(1,2))


hist(data_q$Area,xlab="Area (Millions KM^2)",main="",prob=T,ylim=c(0,.00002),breaks="Scott",axes=FALSE)

density_a<-density(data_q$Area)


polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)


axis(2,at=c(0,.00001,.00002))
axis(1,at=c(0,1000000,3000000,5000000,7000000),labels=c("0","1","3","5","7"))

#abline(v=median(data_q$Area),lty=3,col="red")

#abline(v=mean(data_q$Area),lty=3,col="blue")


hist(data_q$logArea,xlab="ln(Area)",main="",prob=T,ylim=c(0,.17),axes=F)


density_loga<-density(data_q$logArea)


polygon(density_loga,col="black",density=3,angle=45,border="black")

axis(2,at=c(0,.05,.1,.15))


dev.off()


############## Footnote Box-Cox Transform ##################



print("Box Cox Test")

lam<-powerTransform(data_q$Area)
lam<-coef(lam)


### This Gives the result from Footnote TKTK ###

print(lam)




